


#helper functions
computeDistancesM <-function(dist2cent,method,x,mobileCenters,kMov,fixedCenters,kFix,Smatrices){
	#inv <-ginv
	inv <-solve
	for(k in 1:kMov){
		#invCov <- inv(Smatrices[[k]])
		invCov <- (Smatrices[[k]])
		dist2cent[,k] = sqrt(mahalanobis(x,mobileCenters[k,],cov=invCov,inverted=TRUE))
	}
	return(dist2cent)
	
}
computeDistancesM <- cmpfun(computeDistancesM)

computeDistancesE <-function(dist2cent,method,x,mobileCenters,kMov,fixedCenters,kFix,Smatrices){
	for(k in 1:kMov) {
		dist2cent[,k] = sqrt(rowSums(sweep(x,2,mobileCenters[k,],"-")^2))
	}
	return(dist2cent)
}
computeDistancesE <-cmpfun(computeDistancesE)


#One run of vegclust   
vegclustone <-
		function(x,mobileCenters, fixedCenters=NULL, method="NC", m=2,dnoise=NULL, eta = NULL, alpha=0.001, iter.max=100) {
	METHODS <- c( "FCM","NC","GK","GKN")
	method <- match.arg(method, METHODS)
#		if(method=="KM") {
#			m=1.0
#			dnoise=NULL
#			eta=NULL
#		}
	if(method=="FCM" || method=="GK") {
		dnoise=NULL
		eta=NULL
	}
	else if(method=="NC" || method=="GKN") {
		if(is.null(dnoise) || is.na(dnoise)) stop("Must provide a value for dnoise")
		eta = NULL
	}
#		else if(method=="PCM") {
#			if(is.null(eta)) stop("Must provide a vector of values for eta")
#			dnoise = NULL
#		}
	
	x = as.matrix(x)
	n = nrow(x)
	
	#Sets the starting mobile centroids
	#print("Before format choose")
	#print(paste("MCenters",class(mobileCenters),ncol(mobileCenters),sep="::"))
	
	if(is.data.frame(mobileCenters) || is.matrix(mobileCenters)) {
#		print("MCent is data.frame or matrix")
		mobileCenters = as.matrix(mobileCenters)
		if(ncol(mobileCenters)!=ncol(x)) {
			stop("The number and identity of species for mobile centers must be the same as the number and identity of species for x")
		}		
	}
	else if(is.vector(mobileCenters) && length(mobileCenters)==1 && is.numeric(mobileCenters)) {
		if(mobileCenters==1) mobileCenters = t(as.matrix(x[sample(n,mobileCenters),]))
		else mobileCenters = as.matrix(x[sample(n,mobileCenters),])
	}
	else if(is.vector(mobileCenters) && is.numeric(mobileCenters)) {
#		print("Mobile Centers")
#		print(mobileCenters)
#		print(class(mobileCenters))
		mobileCenters = as.matrix(x[mobileCenters,,drop=FALSE])
	} 
	#print("Decision stop")
	if(!is.matrix(mobileCenters)) {
		stop("Provide a number, a vector of seeds, or coordinates for mobile centers")
	}	
	kMov = nrow(mobileCenters)
	rownames(mobileCenters)<-c(1:kMov)
	
	#Sets the fixed centroids
#		if(!is.null(fixedCenters)) {
#			if(is.data.frame(fixedCenters)) {
#				fixedCenters = as.matrix(fixedCenters)
#				kFix = nrow(fixedCenters)
#			}
#			else if(!is.matrix(fixedCenters)) {
#				stop("Fixed centers must be specified as a matrix or a data frame")
#			}	
#			else {kFix = nrow(fixedCenters)}
#			if(ncol(fixedCenters)!=ncol(x)) {
#				stop("The number and identity of species for fixed centers must be the same as the number and identity of species for x")
#			}
#		} else {
#			kFix = 0
#		}
	kFix=0
	NansInSM <-FALSE
	
	#if(method=="PCM" && length(eta)!=(kMov+kFix)) stop("Vector of reference distances (eta) must have a length equal to the number of clusters")
	dist2cent = matrix(0,nrow=n,ncol={kMov+kFix})
	
	if(method=="NC" || method=="GKN") {
		u = matrix(0,nrow=n,ncol={kMov+kFix+1})
		uPrev = matrix(0,nrow=n,ncol={kMov+kFix+1})
	} else {
		u = matrix(0,nrow=n,ncol={kMov+kFix})
		uPrev = matrix(0,nrow=n,ncol={kMov+kFix})
	}
	
	#print("InitSmatrices")
	#initialize S matrices list
	if(method=="GK" || method=="GKN"){
		Smatrices <- lapply(1:{kMov+kFix},FUN=function(n,SmDim){diag(rep(1,SmDim))},ncol(mobileCenters))
	}
	
	continue = TRUE
	iter = 1
	#helper function computes distances
	if(method=="GK" || method == "GKN"){
		computeDistances <- computeDistancesM
	}else{
		computeDistances <- computeDistancesE
	}

	#iterates until no change in memberships or numerical error appears
	while(continue) {
		#		print("ComputeDistances")
		#1. compute distance to centroids (fixed and mobile)
		tmpDist <- computeDistances(dist2cent,method,x,mobileCenters,kMov,fixedCenters,kFix,Smatrices)
		
		if(sum( is.na(tmpDist))>0){
			print(paste("Method: ",method,sep=" "))
			print("NaNs in dist2cent")
			continue=FALSE
		}else{
			dist2cent <- tmpDist
		}
		
		if(continue){
			#2. compute membership to centroids
			if(method=="NC" || method=="GKN") {
				d2cm2<-cbind(dist2cent,dnoise)^2
			}
			else {
				d2cm2<-dist2cent^2
			}
			
			if(sum( is.na(d2cm2))>0){
#				print(paste("Method: ",method,sep=" "))
#				print(paste("distMIN: ",min(dist2cent),sep=" "))
#				print(paste("distMAX: ",max(dist2cent),sep=" "))
#				print(paste("Dnoise NANS?: ",sum(is.na(dnoise)),sep=" "))
#				print(paste("Mcent NANS?: ",sum(is.na(mobileCenters)),sep=" "))
#				print("NaNs in d2cm2")
				continue=FALSE
			}
		}
		if(continue){
			for(k in 1:ncol(d2cm2)) {
				a<-sweep(d2cm2,1,d2cm2[,k],"/")
				u[,k] = 1/rowSums(a^{-1/{m-1}})
			}
			u[d2cm2==0]=1
		}
		#u[is.na(u)] =0
		
		#Check for stopping
		if(iter>2 && continue) {
			continue = {max(abs(u-uPrev))>alpha} && {iter<=iter.max} && {max(abs(u-uPrev2))>alpha}
		}

		if(continue) {
			#3. update mobile centroids
			for(k in 1:kMov) mobileCenters[k,]={{u[,k]^m}%*%x}/sum(u[,k]^m)
			uPrev2 = uPrev
			uPrev = u	   	
			iter=iter+1
			#if(verbose) cat(".")
			if(method =="GK" || method == "GKN"){
				#GK - S computation 
				for(k in 1:{kMov+kFix}){
					tmp1 = sweep(x,2,mobileCenters[k,],"-")
					tmp2 =	sweep(tmp1,1,u[,k]^m,"*")
					Stmp = {t(tmp2) %*% tmp1}/sum(u[,k]^m)
					
					
			#eigs correction
#			StmpRows = nrow(Stmp)
#			tryCatch({
#				ei =eigen(Stmp)
#				if(sum(ei$values[1]/abs(ei$values) >1e15)){
#					ei$values[ei$values[1]/ei$values >1e15] = ei$values[1]/1e15
#					Stmp = ei$vectors %*% diag(ei$values) %*% solve(ei$vectors)
#					
#				}
#			},error=function(e){}) 
				
			DStmp = tryCatch({det(Stmp)},error=function(e){0})
			S<-Smatrices[[k]]
			if(abs(DStmp) > 1e-20){
				#if Stmp can not be inverted do not change Smatrices[[k]]
				tryCatch({S = {DStmp^{1/StmpRows}}*solve(Stmp)},
						error = function(e){
#							print(e$message)
#							print("eigs");print(eigen(Stmp)$values)
#							print("DET");print(det(Stmp))
#							print("Symmetric?");print(isSymmetric(Stmp))
#							#stop("STTOP")
#							warning("Singular matrix")
							
						}
				)
				
			}
			if( sum(is.na( S ) ) >0 ){
				#print("NaNs in Smatrix")
			}else{
				Smatrices[[k]]=S
			}
			
			
				
					
				}
				
			}
			#continue = !NansInSM
		}
		#end of while(continue)
		
	}#end of while
	#if(verbose) cat(paste("\nIterations:", iter,"\n"))
	if(method=="FCM" || method=="KM" || method=="GK") functional = sum((dist2cent^2)*(u^m))
	else if(method=="NC" || method=="GKN") functional = sum({dist2cent^2}*(u[,-{kMov+kFix+1}]^m))+sum(dnoise^2*u[,kMov+kFix+1]^m)
#		else if(method =="GK"){
#			functional = sum( sweep( {dist2cent^2}*{u^m},2,alphas^{m-1},"*" ) )
#		}
#		else if(method =="GKN"){
#			functional = sum( sweep( {dist2cent^2}*{u[,-{kMov+kFix+1}]^m},2,alphas^{m-1},"*" ) )
#							+sum(dnoise^2*u[,kMov+kFix+1]^m)
#		}
	
	
	#Prepare output
	mobileCenters = as.data.frame(mobileCenters)
	u = as.data.frame(u)   
	dist2cent = as.data.frame(dist2cent)   
	for(k in 1:kMov) {
		row.names(mobileCenters)[k] = paste("M",k,sep="")
		names(u)[k] = paste("M",k,sep="")
		names(dist2cent)[k] = paste("M",k,sep="")
	}
#		if(kFix>1) {
#			for(k in {kMov+1}:{kMov+kFix}) {
#				names(u)[k] = paste("F",k,sep="")
#				names(dist2cent)[k] = paste("F",k,sep="")
#			}
#		}
	if(method=="NC" || method=="GKN") names(u)[kMov+kFix+1] = "N"
	rownames(u) = rownames(x)
	rownames(dist2cent) = rownames(x)
	size = colSums(u[,1:{kMov+kFix}])
	withinss = colSums({dist2cent^2}*{u[,1:{kMov+kFix}]^m})
	res = list(mode="raw", method=method, m = m, dnoise = dnoise,eta = eta, memb=u,mobileCenters=mobileCenters, fixedCenters=fixedCenters, dist2clusters=dist2cent, withinss = withinss, size=size, functional=functional, iter=iter)
	class(res)<-"vegclust"
	return(res)
}#end of vegclustone
vegclustone <-cmpfun(vegclustone)



vegclust <-
		function(x,mobileCenters, fixedCenters=NULL, method="NC", m=2,dnoise=NULL, eta = NULL, alpha=0.001, iter.max=100, nstart=1, maxminJ=10, seeds = NULL, verbose=FALSE) {
	
	#print(paste("Vclust start M centers: ",mobileCenters,sep=":::"))
#	if( {sum( is.na(x) )>0} || {sum( is.infinite(x) )>0}  ){
#		stop("Set to score Contains Na or NaN or Inf: in vegclust")
#	}
	if(is.null(seeds)) seeds = 1:nrow(x)
	#print(seeds)
	#If mobileCenters is a number and nstart>1 perform different random starts
	if(is.vector(mobileCenters) && length(mobileCenters)==1 && is.numeric(mobileCenters)) {
		#print("sample")
		bestRun = vegclustone(x,mobileCenters=x[sample(seeds,mobileCenters),,drop=FALSE], fixedCenters, method, m,dnoise, eta, alpha, iter.max)
		if(nstart>1) {
			minJ = 0
			i = 2
			while(i<=nstart) {
				run = vegclustone(x,mobileCenters=x[sample(seeds,mobileCenters),,drop=FALSE], fixedCenters, method, m,dnoise,eta, alpha, iter.max)
				if(run$functional<=bestRun$functional) {
					bestRun = run
					if(max(abs(run$memb-bestRun$memb))<alpha) {
						minJ = minJ+1
					} else {
						minJ = 1
					}
				}
				if(minJ==maxminJ) {
					i=nstart
					if(verbose) print("Maximum number of minimum J reached. Stopping.")
				}
				i=i+1
			}
		}		
		return(bestRun)
	} else { #Perform a single run
		#print(paste("single Run",class(mobileCenters),ncol(mobileCenters),sep=":"))
		return(vegclustone(x,mobileCenters, fixedCenters, method, m,dnoise, eta, alpha, iter.max))
	}
}
vegclust <-cmpfun(vegclust)


